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Abstract - As the penetration of intermittent energy 
sources grows substantially, loads will be required to play an 
increasingly important role in compensating the fast time- 
scale fluctuations in generated power. Recent numerical 
modeling of thermostatically controlled loads (TCLs) has 
demonstrated that such load following is feasible, but analyt- 
ical models that satisfactorily quantify the aggregate power 
consumption of a group of TCLs are desired to enable con- 
troller design. We develop such a model for the aggregate 
power response of a homogeneous population of TCLs to uni- 
form variation of all TCL setpoints. A linearized model of the 
response is derived, and a linear quadratic regulator (LQR) 
has been designed. Using the TCL setpoint as the control 
input, the LQR enables aggregate power to track reference 
signals that exhibit step, ramp and sinusoidal variations. Al- 
though much of the work assumes a homogeneous popula- 
tion of TCLs with deterministic dynamics, we also propose 
a method for probing the dynamics of systems where load 
characteristics are not well known. 

Keywords - Load modeling; load control; renewable 
energy; linear quadratic regulator. 

1 INTRODUCTION 

AS more renewable power generation is added to 
power systems, concerns for grid reliability increase 
due to the intermittency and non-dispatchability associ- 
ated with such sources. Conventional power generators 
have difficulty in manoeuvering to compensate for the 
variability in the power output from renewable sources. 
On the other hand, electrical loads offer the possibility of 
providing the required generation-balancing ancillary ser- 
vices. It is feasible for electrical loads to compensate for 
energy imbalance much more quickly than conventional 
generators, which are often constrained by physical ramp 
rates. 

A population of thermostatically controlled loads 
(TCLs) is well matched to the role of load following. Re- 
search into the behavior of TCLs began with the work of 
HI and IS, who proposed models to capture the hybid dy- 
namics of each thermostat in the population. The aggre- 
gate dynamic response of such loads was investigated by 
in, who derived a coupled ordinary and partial differen- 



tial equation (Fokker-Planck equation) model. The model 
was derived by first assuming a homogeneous group of 
thermostats (all thermostats having the same parameters), 
and then extended using perturbation analysis to obtain the 
model for a non-homogeneous group of thermostats. In 
fS), a discrete-time model of the dynamics of the temper- 
atures of individual thermostats was derived, assuming no 
external random influence. That work was later extended 
by |l6l to introduce random influences and heterogeneity. 

Although the traditional focus has been on direct load 
control methods that directly interrupt power to all loads, 
recent work in jl3| proposed hysteresis-based control by 
manipulating the thermostat setpoint of all loads in the 
population with a common signal. While it is difficult to 
keep track of the temperature and power demands of in- 
dividual loads in the population, the probability of each 
load being in a given state (ON - drawing power or OFF - 
not drawing any power) can be estimated rather accurately. 
System identification techniques were used in |f3]| to obtain 
an aggregate linear TCL model, which was then employed 
in a minimum variance control law to demonstrate the load 
following capability of a population of TCLs. 

In this paper, we derive a transfer function relating the 
aggregate response of a homogeneous group of TCLs to 
disturbances that are applied uniformly to the thermostat 
setpoints of all TCLs. We start from the hybrid temper- 
ature dynamics of individual thermostats in the popula- 
tion, and derive the steady-state probability density func- 
tions of loads being in the ON or OFF states. Using these 
probabilities we calculate aggregate power response to a 
setpoint change. We linearize the response and design a 
linear quadratic regulator to achieve reference tracking by 
the aggregate power demand. WTiile our analytical model 
assumes a homogeneous population of loads, numerical 
studies are proposed to explore situations where there is 
noise and heterogeneity. 

2 STEADY STATE DISTRIBUTION OF LOADS 

2.1 Model development 

The dynamic behavior of the temperature 6{t) of a 
thermostatically controlled cooling-load (TCL), in the ON 
and OFF state and in the absence of noise, can be modeled 



byHa, 
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where 6amb is the ambient temperature, C is the ther- 
mal capacitance, R is the thermal resistance, and P is the 
power drawn by the TCL when in the ON state. This re- 
sponse is shown in Figure [T] 
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Figure 1: Dynamics of temperature of a thermostatic load. 

In steady state the cooling period drives a load from 
temperature 9+ to temperature 9^ . Thus solving ([T]i with 
initial condition 9{Q) = 9^ gives 



9it) = i9a,nb - PR) +9+e~^ . (2) 

From (|2]i we can calculate the steady state cooling time Tc 
by equating 9{Tc) to 9^, 

I PR + 0+ — 9amb 



A similar calculation for the heating time gives, 

Tfe = Cflln(' ^°'"''~^ V (4) 



^amb 



In general, the expressions for the times tc{9f) and tji{9f) 
taken to reach some intermediate temperature 9f during 
the cooling and heating periods, respectively, are. 



t^{9f) ^ CRln 
th{9f) = CRln 
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For a homogeneou^ set of TCL in steady state, the 
number of loads in the ON and OFF states, Nc and Nh 
respectively, will be proportional to their respective cool- 
ing and heating time periods Tc and Th- In the absence 
of any appreciable noise, which ensures that all the loads 
are within the temperature deadband, Nh + Nc — N, we 
obtain. 
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By analogy, it follows that the number of ON-loads 
nc{9) within a temperature band of [9, 9+] is proportional 
to the time taken tc (9) to cool a load down from 9+ to an 
arbitrary temperature 9 > 9^, 
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where ([8]l was used to obtain Likewise, 
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We will denote the ON probability density function by 
and the OFF probability density function by /o(6'), 
while the corresponding cumulative distribution functions 
are denoted Fi{9) and Fo{9), respectively. It is to be 
noted that, Fq (9) is the probability of a load being in OFF 
state and having a temperature 9 E [9- , 9] while Fi (9) 
is the probability of a load being in ON state and having 
a temperature 9 e [9^,9]. Thus, ^0(6*) nh{9)/N and 
Fi{9) = {Nc - nc[9))/N. We can therefore write. 
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and 
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2.2 Simulation 



Figure |2] shows a comparison of the densities calcu- 
lated using ( fTTT ) and ( fT2b and those computed from ac- 
tual simulation of the dynamics of a population of 10,000 
TCLs that included a small amount of noise. The result 
suggests that the assumptions underlying ( fTTT ) and ( fT2] ) are 
reaUstic. 

e ^=32°G, 9 =20 °C, A=2°C 



(actual) 

f^j (actual) 

(analytical) 

L (analytical) 




Figure 2: Steady state densities. 



All loads share the same values for parameters Bamb' C, R and P. 



3 SETPOINT VARIATION 

Control of active power can be achieved by making a 
uniform adjustment to the temperature setpoint of all loads 
within a large population fS). It is assumed that the tem- 
perature deadband moves in unison with the setpoint. Fig- 
ure |3] shows the change in the aggregate power consump- 
tion of a population of TCLs for a small step change in 
the setpoint of all devices. The resulting transient varia- 
tions in the OFF-state and ON-state distributions for the 
population are shown in Figure |4] 
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Figure 3: Change in aggregate power consumption due to a step ctiange 
in temperature setpoint. 

Probability density tunction (OFF state) 



with the deadband width A^9'^ — d'^ =9^ — 9^ remain- 
ing unchanged. The setpoint is shifted by 5 = 6*- — = 
6^ — 9'^. To solve for the power consumption, we need to 
consider four different TCL starting conditions immedi- 
ately after the step change in setpoint, i.e. a-d in Figure|5] 
Using Laplace transforms, we compute the time depen- 
dence of the power consumption for each of these loads 
(shown in Figure|3]l and then compute the total power con- 
sumption by integrating over the distributions /o and fi . 
At the instant the step change in applied, the temperatures 
of loads at points a, b, c and d are 9a, 9b, 9c and 9d, respec- 
tively. 




Figure 5: Different points of interest on the density curves. 
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(a) OFF-state distribution. 
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(b) ON-state distribution. 
Figure 4: Variation in distribution of loads due to setpoint disturbance. 

The aggregate power consumption at any instant in 
time is proportional to the number of loads in the ON state 
at that instant. The first step in quantifying the change in 
power due to a step change in setpoint is therefore to an- 
alyze the behavior of the TCL probability distributions. 
Figure |5] depicts a situation where the setpoint has just 
been increased. The original deadband ranged from 0° 
to 0\, with the setpoint at (0° + 9\)/2. After the posi- 
tive step change, the new deadband lies between 6'_ to 9^, 



ga(t) 
P 













Tc 


Th 










04 {6+} {e.} time 



(a) Power waveform at point a. 
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(c) Power waveform at point c. 
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(d) Power waveform at point d. 
Figure 6: Power waveforms at four different points marked in Figure|5] 

The power consumption ga [t, Ta) of the load at a start- 
ing from the instant when the step change in setpoint is 
applied is shown in Figure [6(a)| All the loads in the OFF- 
state and having a temperature between 6- and B'^ at the 
instant when the deadband shift occurs will have power 
waveforms similar in nature to ga{t, Ta). Thus the load at 
a typifies the behavior of all the loads lying on the OFF- 
state density curve between 6'_ and 6*" . The same argu- 
ment applies for loads at points b, c and d. Figures \6(a)\ 
|6(d)| illustrate the general nature of the power waveforms 
of the loads in all four regions, marked by a, b, c and d in 
Figure [5] 

The Laplace transform of ga {t, Ta) is 



Gais,Ta) 



'G(.) 



where 



G(.) = 



s(l_e-.s(T<,+T,0) 

and Ta = Th - th{Oa), with th{0a) given by Q. Aver- 
aging over all such loads (represented by a) on the OFF 
density curve between temperatures 6'_ and 6'" , we obtain 
the Laplace transform of the average power demand, 



foiOa)Gais,Ta) 



(13) 



where fo{da) can be computed from ( fTTT i. 

In Figure |6(b)[ a load at point b on the ON density 
curve in Figure|5]has power consumption gb{t, Tf,), where 
Th = Tc — tc{9b), and tc{Ob) is given by (|5]i. The Laplace 
transform is. 
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We can compute the average power demand of all the 
loads represented by 6 as 



(14) 



P6(s) = / .fii9b)Gb{s,Tb)d9b 



In Figure [6(c)] the power consumption gc{t,Tc) of a 
load at point c on the OFF density curve in Figure |5] has 
the Laplace transform 



where = Ci? In ( 



The average power de- 
mand of the loads represented by the point c is then given 



by 
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Figure [6(d)] depicts the situation of a load at point d on 
the ON density curve, that suddenly switches to the OFF 
state as the deadband is shifted (for now we assume the 
deadband is shifted to the right, i.e., there is an increase 
in the setpoint). The power consumption gd{t) has the 
Laplace transform 

Gd(s,Td) =e-^(^"+^'')G(s), 

where the dynamics in ([T]i can be solved for = 

CR In I n°'"^''Zg'^ ) ■ The average power demand of the 



loads characterized by point d in Figure |5] is then given 
by 
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(16) 



The average power demand of the whole population 
becomes, 

Pavg{s) = P„(.S) + Pb(.s) + P,(.s) + V d{s) . (17) 

Using (HJI ), ( fT4l i. ( fTSl ) and ( fT6] l we obtain an expression 
for Pat,g(s) that is rather complex. It is hard, and perhaps 
even impossible, to obtain the inverse Laplace transform. 
However, with the assistance of MATHEMATICA®, 
"Pavgis) may be expanded as a series in s. We also make 
use of the assumptions, 

A « [Os - Oarnb + PR) 
A <C {Oarnb — Q s) 
(5 < A 

where 9s is the setpoint temperature. Note that the first 
two assumptions require that the deadband width is small, 
while the third assumption requires that the shift in the 
deadband is small relative to the deadband width. This 
latter assumption ensures that the load densities are not 
perturbed far from their steady-state forms. Accordingly, 
the steady-state power consumption is given by 

p _ [Oamb - e+)N 

ijR 

where rj is the electrical efficiency of the cooling equip- 
ment and is the population size. The deviation in power 
response can be approximated by 



s s 



(18) 



where 
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and Tco and Tho are the original (prior to the setpoint shift) 
steady-state cooUng and heating times, respectively, given 
by (O and (HJi. The transfer function for this linear model 

is. 



T{s) 



S/s 



= - d 



Due to the assumptions of low-noise and homogeneity, 
our analytical model is undamped. The actual system, on 
the other hand, experiences both heterogeneity and noise, 
and therefore will exhibit a damped response. In order to 
capture that effect, we have chosen to add a damping term 
<T (to be estimated on-line) into the model, giving 
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(19) 



Figure |3] shows a comparison between the response 
calculated from the model ( fT9] l and the true response to 
a step change in the setpoint obtained from simulation. A 
damping coefficient of 0.002 min^^ was added, as that 
value gave a close match to the decay in the actual system 
response. 



Total Power Variation 




where yd is the reference, as the third state w{t) = 
/o (?/(''') ~ yd{T))dT of the system, the modified state- 
space model becomes 



i = Ax + Bu + E yd 
y = C_x_ + Dm 



where x =\x and. 
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Minimizing the cost function 



J= I {x{t)^Qx{t)+u{tfR)dt 
'o 



Figure 7: Comparison of the approximate model witli the actual simula- 
tion, for the same setpoint disturbance as in Figure|3] 



4 CONTROL LAW 

The TCL load controller, described by the transfer 
function (fT9] l, can also be expressed in state-space form. 
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where the input u{t) is the shift in the deadband of all 
TCLs, and the output y{t) is the change in the total power 
demand from the steady-state value. The state-space ma- 
trices are given by 
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C = [-1 0] , 
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Our goal is to design a controller using the linear 
quadratic regulator (LQR) approach Q to track an ex- 
ogenous reference yd- We observe that the system has 
an open-loop zero very close to the imaginary axis {d ^ 
ojAa) and hence we need to use an integral controller 
Considering the integral of the output error e = [y — yd). 



where Q > Osxs ™d i? > are design variables, we 
obtain the optimal control law u{t) of the form 



u = -(Kx + Gyd), 



with G a pre-compensator gain chosen to ensure unity DC 
gain. Since we can only measure the output y{t) and the 
third state w{t), the other two states are estimated using 
a linear quadratic estimator Q which has the state-space 
form. 



x = Ax + Bw + 'L{y — yd) 
ij = Cx + Dm 



u = -K 



Gyd- 



Total Power Variation Temperature Set-point Variation Figure 8: Reference tracking achieved througli setpoint sliift 




(a) Response to step reference and tlie control input 



(a) OFF-state distribution. 
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(b) Response to ramp reference and tlie control input 
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(b) ON-state distribution. 



(c) Response to sinusoidal reference and the control input 



Figure 9: Variation in distribution of loads under the influence of the 
controller. 



The plots in Figure [8] show that the controller can be 
used to force the aggregate power demand of the TCL pop- 
ulation to track a range of reference signals. The transient 
variations in the ON-state and OFF-state populations are 
shown in Figure |9] In comparison with the uncontrolled 
response of Figure IH it can be seen that the controller 
suppresses the lengthy oscillations. Figure |9] shows that 
in presence of the controller, the distribution of loads al- 
most always remains close to steady state, justifying an 
assumption made during the derivation of the model. 



5 HETEROGENEITY AND NOISE 



The work presented in previous sections assumes a ho- 
mogeneous population of loads with deterministic dynam- 
ics. The analysis remains valid if we consider the pos- 
sibility of grouping a large number loads having closely 
matched parameters, and if there is very low noise in the 
system. When such assumptions no longer remain valid, 
we cannot design a tracking controller based on the devel- 
oped model. In such cases, however, we propose a probing 
method that can balance over- or under-production of en- 



ergy over a certain duration of time. 
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(a) Power response to short duration pulses. 
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(b) Energy consumed (over nominal) in response to such pulses. (A neg- 
ative value means energy is delivered.) 

Figure 10: Energy consumption in the probing method. 

Figure |10(a)| can be used to explain the probing 
method. The temperature setpoint is increased and held at 
that value for a short duration and then returned to its orig- 
inal value. The system is probed by short pulses spaced 
reasonably far from each other in time. The energy deliv- 
ered during such probing is monitored, with Figure [TO(b)| 
providing an illustration. It can be seen that the energy 
consumed, relative to the nominal consumption, is actu- 
ally negative suggesting that energy is "delivered" by the 
loads when probed with positive pulses. Knowing that 
over a certain duration a certain amount of energy can be 
delivered by the loads, the pulses can be scheduled to bal- 
ance any under-generation. Similarly, over-generation can 
be balanced using negative pulses. 

6 CONCLUSION 

In this paper we have analytically derived a transfer 
function relating the change in aggregate power demand 
of a population of TCLs to a change in thermostat setpoint 
applied to all TCLs in unison. We have designed a linear 
quadratic regulator to enable the aggregate power demand 



to track reference signals. This suggests the derived aggre- 
gate response model could be used to allow load to track 
fluctuations in renewable generation. The analysis has 
been based on the assumptions that the TCL population 
is homogeneous and that the noise level is insignificant. 
When such assumptions do not hold, we propose a prob- 
ing method that can be used to perform energy balance. 
Further studies are required to incorporate the effects of 
heterogeneity and noise into the model. Those extensions 
are important for determining the damping coefficient. 

Similar analysis can be used to establish the aggregate 
characteristics of groups of plug-in electric vehicles, an- 
other candidate for compensating the variability in renew- 
able generation. 
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